knitr::opts_chunk$set(echo = TRUE)

cmdstanr::set_cmdstan_path(path = "C:/Users/kueng/.cmdstan/cmdstan-2.35.0")

library(tidyverse)
library(R.utils)
library(wbCorr)
library(readxl)
library(kableExtra)
library(brms)
library(bayesplot)
library(see)
library(beepr)
library(DHARMa)
library(digest)



source(file.path('Functions', 'ReportModels.R'))
source(file.path('Functions', 'PrettyTables.R'))
source(file.path('Functions', 'ReportMeasures.R'))
source(file.path('Functions', 'PrepareData.R'))

report_function_hash <- digest::digest(summarize_brms)
system("shutdown /a")
## [1] 1116
# Set options for analysis
use_mi = FALSE
shutdown = FALSE
report_ordinal = FALSE
report_hurdle = TRUE
do_priorsense = FALSE

options(
  dplyr.print_max = 100, 
  brms.backend = 'cmdstan',
  brms.file_refit = ifelse(use_mi, 'never', 'on_change'),
  brms.file_refit = 'on_change',
  #brms.file_refit = 'always',
  error = function() {
    beepr::beep(sound = 5)
    if (shutdown) {
      system("shutdown /s /t 180")
      quit(save = "no", status = 1)
    }
  }
  , es.use_symbols = TRUE
)


####################### Model parameters #######################

iterations = 12000 # 10'000 per chain to achieve 40'000
warmup = 2000

# NO AR!!!
#corstr = 'ar'
#corstr = 'cosy_couple'
#corstr = 'cosy_couple:user'


################################################################

suffix = as.character(iterations)
df <- openxlsx::read.xlsx(file.path('long.xlsx'))
df_original <- df

df_double <- prepare_data(df, recode_pushing = TRUE, use_mi = use_mi)[[1]]

Constructing scales Re-coding pusing reshaping data (4field) centering data within and between

summary(df_double$pushing)

Min. 1st Qu. Median Mean 3rd Qu. Max. NA’s 0.0000 0.0000 0.0000 0.1635 0.0000 5.0000 241

Modelling

# For indistinguishable Dyads
model_rows_fixed <- c(
    'Intercept', 
    # '-- WITHIN PERSON MAIN EFFECTS --', 
    'persuasion_self_cw', 
    'persuasion_partner_cw', 
    'pressure_self_cw', 
    'pressure_partner_cw', 
    'pushing_self_cw', 
    'pushing_partner_cw', 
    'day', 
    'planPlan',
    'weartime_self_cw',
    
    # '-- BETWEEN PERSON MAIN EFFECTS',
    'persuasion_self_cb',
    'persuasion_partner_cb',
    'pressure_self_cb',
    'pressure_partner_cb',
    'pushing_self_cb',
    'pushing_partner_cb',
    'weartime_self_cb'
  )


model_rows_fixed_ordinal <- c(
  model_rows_fixed[1],
  'Intercept[1]',
  'Intercept[2]',
  'Intercept[3]',
  'Intercept[4]',
  'Intercept[5]',
  model_rows_fixed[2:length(model_rows_fixed)]
)

model_rows_random <- c(
  # '--------------',
  # '-- RANDOM EFFECTS --',
  'sd(Intercept)', 
  'sd(persuasion_self_cw)',
  'sd(persuasion_partner_cw)',
  'sd(pressure_self_cw)',
  'sd(pressure_partner_cw)',
  'sd(pushing_self_cw)',
  'sd(pushing_partner_cw)',
  # '-- CORRELATION STRUCTURE -- ', 
  'sigma'
)

model_rows_random_ordinal <- c(model_rows_random,'disc')
# For indistinguishable Dyads
model_rownames_fixed <- c(
    "Intercept", 
    # "-- WITHIN PERSON MAIN EFFECTS --", 
    "Daily persuasion experienced", 
    "Daily persuasion utilized (partner's view)", # OR partner received
    "Daily pressure experienced", 
    "Daily pressure utilized (partner's view)", 
    "Daily pushing experienced", 
    "Daily pushing utilized (partner's view)", 
    "Day", 
    "Actionplan",
    "Daily weartime",
    
    # "-- BETWEEN PERSON MAIN EFFECTS",
    "Mean persuasion experienced", 
    "Mean persuasion utilized (partner's view)", 
    "Mean pressure experienced", 
    "Mean pressure utilized (partner's view)", 
    "Mean pushing experienced", 
    "Mean pushing utilized (partner's view)", 
    "Mean weartime"
  )


model_rownames_fixed_ordinal <- c(
  model_rownames_fixed[1],
  'Intercept[1]',
  'Intercept[2]',
  'Intercept[3]',
  'Intercept[4]',
  'Intercept[5]',
  model_rownames_fixed[2:length(model_rownames_fixed)]
)

model_rownames_random <- c(
  # '--------------',
  # '-- RANDOM EFFECTS --',
  'sd(Intercept)', 
  "sd(Daily persuasion experienced)", 
  "sd(Daily persuasion utilized (partner's view))", # OR partner received
  "sd(Daily pressure experienced)", 
  "sd(Daily pressure utilized (partner's view))", 
  "sd(Daily pushing experienced)", 
  "sd(Daily pushing utilized (partner's view))", 
  # '-- CORRELATION STRUCTURE -- ', 
  'sigma'
)

model_rownames_random_ordinal <- c(model_rownames_random,'disc')
rows_to_pack <- list(
  "Within-Person Effects" = c(2,10),
  "Between-Person Effects" = c(11,17),
  "Random Effects" = c(18, 24), 
  "Additional Parameters" = c(25,25)
  )


rows_to_pack_ordinal <- list(
  "Intercepts" = c(1,6),
  "Within-Person Effects" = c(2+5,10+5),
  "Between-Person Effects" = c(11+5,17+5),
  "Random Effects" = c(18+5, 24+5), 
  "Additional Parameters" = c(25+5,25+6)
  )

HURDLE MODELS

# For indistinguishable Dyads
model_rows_fixed_hu <- c(
    'Intercept', 
    'hu_Intercept',
    # '-- WITHIN PERSON MAIN EFFECTS --', 
    'persuasion_self_cw', 
    'persuasion_partner_cw', 
    'pressure_self_cw', 
    'pressure_partner_cw', 
    'pushing_self_cw', 
    'pushing_partner_cw', 
    'day', 
    'planPlan',
    'weartime_self_cw',
    
    # '-- BETWEEN PERSON MAIN EFFECTS',
    'persuasion_self_cb',
    'persuasion_partner_cb',
    'pressure_self_cb',
    'pressure_partner_cb',
    'pushing_self_cb',
    'pushing_partner_cb',
    'weartime_self_cb',
    
    # HURDLE MODEL
    # '-- WITHIN PERSON MAIN EFFECTS --', 
    'hu_persuasion_self_cw', 
    'hu_persuasion_partner_cw', 
    'hu_pressure_self_cw', 
    'hu_pressure_partner_cw', 
    'hu_pushing_self_cw', 
    'hu_pushing_partner_cw', 
    'hu_day', 
    'hu_planPlan',
    'hu_weartime_self_cw',
    
    # '-- BETWEEN PERSON MAIN EFFECTS',
    'hu_persuasion_self_cb',
    'hu_persuasion_partner_cb',
    'hu_pressure_self_cb',
    'hu_pressure_partner_cb',
    'hu_pushing_self_cb',
    'hu_pushing_partner_cb',
    'hu_weartime_self_cb'
  )

model_rows_fixed_hu_ordinal <- c(
  model_rows_fixed_hu[1:2],
  'Intercept[1]',
  'Intercept[2]',
  'Intercept[3]',
  'Intercept[4]',
  'Intercept[5]',
  model_rows_fixed_hu[3:length(model_rows_fixed_hu)]
)


model_rows_random_hu <- c(
  # '--------------',
  # '-- RANDOM EFFECTS --',
  'sd(Intercept)', 
  'sd(hu_Intercept)',
  'sd(persuasion_self_cw)',
  'sd(persuasion_partner_cw)',
  'sd(pressure_self_cw)',
  'sd(pressure_partner_cw)',
  'sd(pushing_self_cw)',
  'sd(pushing_partner_cw)',
  # HURDLE
  'sd(hu_persuasion_self_cw)',
  'sd(hu_persuasion_partner_cw)',
  'sd(hu_pressure_self_cw)',
  'sd(hu_pressure_partner_cw)',
  'sd(hu_pushing_self_cw)',
  'sd(hu_pushing_partner_cw)',
  # '-- CORRELATION STRUCTURE -- ', 
  'sigma'
)

model_rows_random_hu_ordinal <- c(model_rows_random_hu,'disc')
# For indistinguishable Dyads
model_rownames_fixed_hu <- c(
    "Intercept", 
    "Hurdle Intercept",
    # "-- WITHIN PERSON MAIN EFFECTS --", 
    "Daily persuasion experienced", 
    "Daily persuasion utilized (partner's view)", # OR partner received
    "Daily pressure experienced", 
    "Daily pressure utilized (partner's view)", 
    "Daily pushing experienced", 
    "Daily pushing utilized (partner's view)", 
    "Day", 
    "Actionplan",
    "Daily weartime",
    
    # "-- BETWEEN PERSON MAIN EFFECTS",
    "Mean persuasion experienced", 
    "Mean persuasion utilized (partner's view)", 
    "Mean pressure experienced", 
    "Mean pressure utilized (partner's view)", 
    "Mean pushing experienced", 
    "Mean pushing utilized (partner's view)", 
    "Mean weartime",
    
    # HURDLE
    # "-- WITHIN PERSON MAIN EFFECTS --", 
    "Hu Daily persuasion experienced", 
    "Hu Daily persuasion utilized (partner's view)", # OR partner received
    "Hu Daily pressure experienced", 
    "Hu Daily pressure utilized (partner's view)", 
    "Hu Daily pushing experienced", 
    "Hu Daily pushing utilized (partner's view)", 
    "Hu Day", 
    "Hu Actionplan",
    "Hu Daily weartime",
    
    # "-- BETWEEN PERSON MAIN EFFECTS",
    "Hu Mean persuasion experienced", 
    "Hu Mean persuasion utilized (partner's view)", 
    "Hu Mean pressure experienced", 
    "Hu Mean pressure utilized (partner's view)", 
    "Hu Mean pushing experienced", 
    "Hu Mean pushing utilized (partner's view)", 
    "Hu Mean weartime"
  )



model_rownames_fixed_hu_ordinal <- c(
  model_rownames_fixed_hu[1:2],
  'Intercept[1]',
  'Intercept[2]',
  'Intercept[3]',
  'Intercept[4]',
  'Intercept[5]',
  model_rownames_fixed_hu[3:length(model_rownames_fixed_hu)]
)



model_rownames_random_hu <- c(
  # '--------------',
  # '-- RANDOM EFFECTS --',
  'sd(Intercept)', 
  'sd(Hurdle Intercept)', 
  "sd(Daily persuasion experienced)", 
  "sd(Daily persuasion utilized (partner's view))", # OR partner received
  "sd(Daily pressure experienced)", 
  "sd(Daily pressure utilized (partner's view))", 
  "sd(Daily pushing experienced)", 
  "sd(Daily pushing utilized (partner's view))", 
  
  # Hurdle
  "sd(Hu Daily persuasion experienced)", 
  "sd(Hu Daily persuasion utilized (partner's view))", # OR partner received
  "sd(Hu Daily pressure experienced)", 
  "sd(Hu Daily pressure utilized (partner's view))", 
  "sd(Hu Daily pushing experienced)", 
  "sd(Hu Daily pushing utilized (partner's view))", 
  # '-- CORRELATION STRUCTURE -- ', 
  'sigma'
)

model_rownames_random_hu_ordinal <- c(model_rownames_random_hu,'disc')
rows_to_pack_hu <- list(
  "Conditional Within-Person Effects" = c(3,11),
  "Conditional Between-Person Effects" = c(12,18),
  
  "Hurdle Within-Person Effects" = c(19,27),
  "Hurdle Between-Person Effects" = c(28,34),
  
  "Random Effects" = c(35, 48), 
  "Additional Parameters" = c(49,49)
  )

rows_to_pack_hu_ordinal <- list(
  "Intercepts" = c(1,7),
  "Conditional Within-Person Effects" = c(3+5,11+5),
  "Conditional Between-Person Effects" = c(12+5,18+5),
  
  "Hurdle Within-Person Effects" = c(19+5,27+5),
  "Hurdle Between-Person Effects" = c(28+5,34+5),
  
  "Random Effects" = c(35+5, 48+5), 
  "Additional Parameters" = c(49+5,49+6)
  )

Self-Reported MVPA

range(df_double$pa_sub, na.rm = T) 
## [1]   0 720
hist(df_double$pa_sub, breaks = 40) 

hist(log(df_double$pa_sub+00000000001), breaks = 40)

Hurdle Lognormal Model

formula <- bf(
  pa_sub ~ 
    persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw +
    
    persuasion_self_cb + persuasion_partner_cb +
    pressure_self_cb + pressure_partner_cb +
    pushing_self_cb + pushing_partner_cb +
    
    plan +
    day + 
    
    # Random effects
    (1 + persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw | coupleID),
  
  hu = ~ persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw +
    
    persuasion_self_cb + persuasion_partner_cb +
    pressure_self_cb + pressure_partner_cb +
    pushing_self_cb + pushing_partner_cb +
    
    plan +
    day + 
    
    # Random effects
    (1 + persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw | coupleID)
  #, autocor = autocor_str
) 

prior1 <- c(
  brms::set_prior("normal(0, 2.5)", class = "b")
  , brms::set_prior("normal(0, 2)", class = "b", dpar = "hu")
  , brms::set_prior("normal(0, 50)", class = "Intercept") # for non-zero PA
  , brms::set_prior("normal(0.5, 2.5)", class = "Intercept", dpar = 'hu') # hurdle part
  , brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
  , brms::set_prior("student_t(3, 0, 2.5)", class = "sigma", lb = 0)
)

#brms::validate_prior(
#  prior1, 
#  formula = formula, 
#  data = df_double, 
#  family = hurdle_lognormal()
#)

#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

pa_sub <- my_brm(
  mi = use_mi, 
  imputed_data = implist,
  
  formula = formula, 
  prior = prior1,
  data = df_double, 
  family = brms::hurdle_lognormal(), 
  #family = brms::hurdle_negbinomial(), 
  #family = brms::hurdle_poisson(),
  #control = list(adapt_delta = 0.95, max_treedepth = 15),
  iter = iterations,
  warmup = warmup,
  chains = 4,
  cores = 4,
  seed = 42,
  file = file.path("models_cache_brms", paste0("pa_sub_hu_lognormal_NOAR", suffix))
  #, file_refit = 'always'
)
## Warning: Rows containing NAs were excluded from the model.
## Start sampling
pa_sub_digest <- digest::digest(pa_sub)
check_brms(pa_sub, log_pp_check = TRUE)
## 
## Divergences:
## 0 of 40000 iterations ended with a divergence.
## 
## Tree depth:
## 0 of 40000 iterations saturated the maximum tree depth of 10.
## 
## Energy:
## E-BFMI indicated no pathological behavior.

## Using 10 posterior draws for ppc type 'ecdf_overlay' by default.

## Using 10 posterior draws for ppc type 'dens_overlay' by default.

## Warning: Found 1 observations with a pareto_k > 0.7 in model 'model'. We
## recommend to set 'moment_match = TRUE' in order to perform moment matching for
## problematic observations.
## 
## Computed from 40000 by 3736 log-likelihood matrix.
## 
##          Estimate    SE
## elpd_loo -10130.3 165.9
## p_loo       175.4   6.1
## looic     20260.6 331.9
## ------
## MCSE of elpd_loo is NA.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.5, 1.9]).
## 
## Pareto k diagnostic values:
##                          Count Pct.    Min. ESS
## (-Inf, 0.7]   (good)     3735  100.0%  2953    
##    (0.7, 1]   (bad)         1    0.0%  <NA>    
##    (1, Inf)   (very bad)    0    0.0%  <NA>    
## See help('pareto-k-diagnostic') for details.
DHARMa.check_brms.all(pa_sub, integer = TRUE, outliers_type = 'bootstrap')

## 
##  DHARMa bootstrapped outlier test
## 
## data:  model.check
## outliers at both margin(s) = 9, observations = 3736, p-value = 0.32
## alternative hypothesis: two.sided
##  percent confidence interval:
##  0.0005353319 0.0032119914
## sample estimates:
## outlier frequency (expected: 0.00159796573875803 ) 
##                                        0.002408994
if (do_priorsense) {
  priorsense_vars <- c(
      'Intercept',
      'b_persuasion_self_cw',
      'b_persuasion_partner_cw',
      'b_pressure_self_cw',
      'b_pressure_partner_cw',
      'b_pushing_self_cw',
      'b_pushing_partner_cw'
  )
  
  hurdle_priorsense_vars <- c(
    'Intercept_hu',
    'b_hu_persuasion_self_cw',
    'b_hu_persuasion_partner_cw',
    'b_hu_pressure_self_cw',
    'b_hu_pressure_partner_cw',
    'b_hu_pushing_self_cw',
    'b_hu_pushing_partner_cw'
  )
  
  gc()
  priorsense::powerscale_sensitivity(pa_sub, variable = c(priorsense_vars, hurdle_priorsense_vars))
  priorsense::powerscale_plot_dens(pa_sub, variable = c(priorsense_vars, hurdle_priorsense_vars))
  priorsense::powerscale_plot_ecdf(pa_sub, variable = c(priorsense_vars, hurdle_priorsense_vars))
  priorsense::powerscale_plot_quantities(pa_sub, variable = c(priorsense_vars, hurdle_priorsense_vars))
}
# summarize with rope range for hurdle part
summary_pa_sub_hurdle <- summarize_brms(
  pa_sub, 
  stats_to_report = c('pd','ROPE'),
  rope_range = c(-0.18, 0.18),
  model_rows_fixed = model_rows_fixed_hu,
  model_rows_random = model_rows_random_hu,
  model_rownames_fixed = model_rownames_fixed_hu,
  model_rownames_random = model_rownames_random_hu,
  exponentiate = T) 
## Warning in summarize_brms(pa_sub, stats_to_report = c("pd", "ROPE"), rope_range
## = c(-0.18, : Coefficients were exponentiated. Double check if this was
## intended.
# rope range for continuous part of the model
rope_factor <- sd(log(pa_sub$data$pa_sub[pa_sub$data$pa_sub > 0]))
rope_range_continuous = c(-0.1 * rope_factor, 0.1 * rope_factor)

summary_pa_sub_continuous <- summarize_brms(
  pa_sub, 
  rope_range = rope_range_continuous,
  model_rows_fixed = model_rows_fixed_hu,
  model_rows_random = model_rows_random_hu,
  model_rownames_fixed = model_rownames_fixed_hu,
  model_rownames_random = model_rownames_random_hu,
  exponentiate = T) 
## Sampling priors, please wait...
## Warning in summarize_brms(pa_sub, rope_range = rope_range_continuous,
## model_rows_fixed = model_rows_fixed_hu, : Coefficients were exponentiated.
## Double check if this was intended.
# Replace only the ROPE and % in Rope columns for rows with 'Hu'
summary_pa_sub <- summary_pa_sub_continuous

columns_to_replace <- c("ROPE", "inside ROPE")

summary_pa_sub[grepl('Hu', rownames(summary_pa_sub)), columns_to_replace] <- 
  summary_pa_sub_hurdle[grepl('Hu', rownames(summary_pa_sub_hurdle)), columns_to_replace]

# Print the updated dataframe
summary_pa_sub %>%
  print_df(rows_to_pack = rows_to_pack_hu)
exp(Est.) SE 95% CI pd ROPE inside ROPE BF BF_Evidence Rhat Bulk_ESS Tail_ESS
Intercept 38.88*** 2.67 [33.95, 44.46] 1.000 [0.92, 1.08] 0.000 >100 Overwhelming Evidence 1.000 9894 20019
Hurdle Intercept 0.27*** 0.04 [ 0.20, 0.36] 1.000 [0.84, 1.20] 0.000 >100 Overwhelming Evidence 1.001 8325 15634
Conditional Within-Person Effects
Daily persuasion experienced 1.03 0.03 [ 0.98, 1.09] 0.862 [0.92, 1.08] 0.965 0.019 Very Strong Evidence for Null 1.000 14778 22659
Daily persuasion utilized (partner’s view) 1.03 0.02 [ 0.99, 1.08] 0.922 [0.92, 1.08] 0.975 0.026 Very Strong Evidence for Null 1.000 18333 24464
Daily pressure experienced 0.91 0.04 [ 0.82, 1.00] 0.972 [0.92, 1.08] 0.356 0.140 Moderate Evidence for Null 1.000 31209 25873
Daily pressure utilized (partner’s view) 0.94 0.04 [ 0.86, 1.03] 0.917 [0.92, 1.08] 0.651 0.046 Strong Evidence for Null 1.000 31187 27274
Daily pushing experienced 0.99 0.03 [ 0.93, 1.06] 0.614 [0.92, 1.08] 0.971 0.013 Very Strong Evidence for Null 1.000 24616 26607
Daily pushing utilized (partner’s view) 0.97 0.03 [ 0.91, 1.03] 0.859 [0.92, 1.08] 0.931 0.021 Very Strong Evidence for Null 1.000 28089 29033
Day 0.99 0.06 [ 0.88, 1.12] 0.559 [0.92, 1.08] 0.794 0.025 Very Strong Evidence for Null 1.000 40169 30324
Actionplan 1.36*** 0.06 [ 1.25, 1.49] 1.000 [0.92, 1.08] 0.000 >100 Overwhelming Evidence 1.000 38568 29578
Daily weartime NA NA NA NA NA NA NA NA NA NA NA
Conditional Between-Person Effects
Mean persuasion experienced 1.04 0.16 [ 0.76, 1.40] 0.588 [0.92, 1.08] 0.380 0.064 Strong Evidence for Null 1.000 6820 13057
Mean persuasion utilized (partner’s view) 0.98 0.15 [ 0.73, 1.33] 0.543 [0.92, 1.08] 0.393 0.061 Strong Evidence for Null 1.000 7023 13637
Mean pressure experienced 1.18 0.21 [ 0.84, 1.69] 0.831 [0.92, 1.08] 0.225 0.112 Moderate Evidence for Null 1.000 10304 19317
Mean pressure utilized (partner’s view) 0.92 0.17 [ 0.64, 1.32] 0.678 [0.92, 1.08] 0.297 0.082 Strong Evidence for Null 1.000 10797 19928
Mean pushing experienced 1.21 0.27 [ 0.78, 1.91] 0.800 [0.92, 1.08] 0.192 0.128 Moderate Evidence for Null 1.001 9138 16846
Mean pushing utilized (partner’s view) 1.33 0.31 [ 0.84, 2.11] 0.893 [0.92, 1.08] 0.124 0.199 Moderate Evidence for Null 1.001 9645 17126
Mean weartime NA NA NA NA NA NA NA NA NA NA NA
Hurdle Within-Person Effects
Hu Daily persuasion experienced 1.59*** 0.11 [ 1.39, 1.85] 1.000 [0.84, 1.20] 0.000 >100 Overwhelming Evidence 1.000 19808 23760
Hu Daily persuasion utilized (partner’s view) 1.34*** 0.09 [ 1.19, 1.54] 1.000 [0.84, 1.20] 0.031 >100 Overwhelming Evidence 1.000 24801 23816
Hu Daily pressure experienced 0.96 0.15 [ 0.70, 1.31] 0.595 [0.84, 1.20] 0.747 0.077 Strong Evidence for Null 1.000 28143 25585
Hu Daily pressure utilized (partner’s view) 1.48* 0.28 [ 1.04, 2.29] 0.987 [0.84, 1.20] 0.122 1.070 Weak Evidence 1.000 27271 19308
Hu Daily pushing experienced 0.95 0.14 [ 0.71, 1.31] 0.630 [0.84, 1.20] 0.744 0.079 Strong Evidence for Null 1.000 21206 24359
Hu Daily pushing utilized (partner’s view) 1.32* 0.15 [ 1.07, 1.68] 0.995 [0.84, 1.20] 0.177 1.912 Weak Evidence 1.000 32065 26478
Hu Day 0.86 0.13 [ 0.64, 1.15] 0.844 [0.84, 1.20] 0.573 0.122 Moderate Evidence for Null 1.000 44632 29875
Hu Actionplan 10.08*** 0.96 [ 8.39, 12.18] 1.000 [0.84, 1.20] 0.000 >100 Overwhelming Evidence 1.000 43764 30032
Hu Daily weartime NA NA NA NA NA NA NA NA NA NA NA
Hurdle Between-Person Effects
Hu Mean persuasion experienced 1.39 0.50 [ 0.68, 2.81] 0.816 [0.84, 1.20] 0.260 0.278 Moderate Evidence for Null 1.001 6046 12697
Hu Mean persuasion utilized (partner’s view) 1.31 0.47 [ 0.64, 2.66] 0.770 [0.84, 1.20] 0.290 0.236 Moderate Evidence for Null 1.001 6006 12353
Hu Mean pressure experienced 0.41* 0.17 [ 0.18, 0.92] 0.985 [0.84, 1.20] 0.036 2.207 Weak Evidence 1.000 10554 19970
Hu Mean pressure utilized (partner’s view) 0.50 0.21 [ 0.22, 1.14] 0.951 [0.84, 1.20] 0.090 0.818 Weak Evidence for Null 1.001 11656 21157
Hu Mean pushing experienced 1.17 0.62 [ 0.41, 3.40] 0.619 [0.84, 1.20] 0.253 0.277 Moderate Evidence for Null 1.000 8908 17485
Hu Mean pushing utilized (partner’s view) 2.19 1.16 [ 0.76, 6.17] 0.928 [0.84, 1.20] 0.093 0.788 Weak Evidence for Null 1.000 8926 17467
Hu Mean weartime NA NA NA NA NA NA NA NA NA NA NA
Random Effects
sd(Intercept) 0.31 0.04 [0.24, 0.41] NA NA NA NA NA 1.000 9257 16013
sd(Hurdle Intercept) 0.74 0.11 [0.57, 1.00] NA NA NA NA NA 1.000 9776 19593
sd(Daily persuasion experienced) 0.11 0.02 [0.08, 0.17] NA NA NA NA NA 1.000 16713 24153
sd(Daily persuasion utilized (partner’s view)) 0.08 0.02 [0.04, 0.13] NA NA NA NA NA 1.000 20429 21942
sd(Daily pressure experienced) 0.07 0.06 [0.00, 0.22] NA NA NA NA NA 1.000 14773 17080
sd(Daily pressure utilized (partner’s view)) 0.06 0.05 [0.00, 0.18] NA NA NA NA NA 1.000 17040 15424
sd(Daily pushing experienced) 0.09 0.04 [0.02, 0.17] NA NA NA NA NA 1.000 11995 9340
sd(Daily pushing utilized (partner’s view)) 0.08 0.03 [0.01, 0.15] NA NA NA NA NA 1.000 11491 8937
sd(Hu Daily persuasion experienced) 0.22 0.09 [0.04, 0.41] NA NA NA NA NA 1.000 8029 7307
sd(Hu Daily persuasion utilized (partner’s view)) 0.18 0.09 [0.02, 0.36] NA NA NA NA NA 1.000 8630 7245
sd(Hu Daily pressure experienced) 0.18 0.17 [0.01, 0.70] NA NA NA NA NA 1.000 14702 16564
sd(Hu Daily pressure utilized (partner’s view)) 0.23 0.21 [0.01, 0.90] NA NA NA NA NA 1.000 13874 18127
sd(Hu Daily pushing experienced) 0.58 0.17 [0.31, 0.99] NA NA NA NA NA 1.000 16458 24605
sd(Hu Daily pushing utilized (partner’s view)) 0.22 0.14 [0.01, 0.54] NA NA NA NA NA 1.000 13039 12440
Additional Parameters
sigma 0.68 0.01 [0.65, 0.70] NA NA NA NA NA 1.000 45378 29070 
# Plot continuous part of model

variable <- c(
  '(Intercept)',
  'b_persuasion_self_cw',
  'b_persuasion_partner_cw',
  'b_pressure_self_cw',
  'b_pressure_partner_cw',
  'b_pushing_self_cw',
  'b_pushing_partner_cw'
)


plot(
  bayestestR::p_direction(pa_sub, parameter = variable),
  priors = TRUE
) + theme_bw()
## Warning in `==.default`(dens$Parameter, parameter): longer object length is not
## a multiple of shorter object length
## Warning in is.na(e1) | is.na(e2): longer object length is not a multiple of
## shorter object length

plot(
  bayestestR::rope(
    pa_sub, 
    parameter = variable, 
    range = rope_range_continuous,
    verbose = F,
    ci = 1
  )
) + theme_bw()

# Hurdle part of the model
variable <- c(
  'b_hu_persuasion_self_cw',
  'b_hu_persuasion_partner_cw',
  'b_hu_pressure_self_cw',
  'b_hu_pressure_partner_cw',
  'b_hu_pushing_self_cw',
  'b_hu_pushing_partner_cw'
)

plot(
  bayestestR::p_direction(pa_sub, parameter = variable),
  priors = TRUE
) + theme_bw()

# The rope range for the bernoulli part of the model is -0.18, 0.18
plot(
  bayestestR::rope(pa_sub, parameter = variable, range = c(-0.18, 0.18), ci = 1),
  verbose = FALSE
) + theme_bw()
## Possible multicollinearity between b_persuasion_partner_cb and
##   b_persuasion_self_cb (r = 0.76), b_hu_persuasion_partner_cb and
##   b_hu_persuasion_self_cb (r = 0.75). This might lead to inappropriate
##   results. See 'Details' in '?rope'.

conds_eff <- conditional_spaghetti(
  pa_sub, 
  effects = c(
    'persuasion_self_cw',
    'persuasion_partner_cw',
    'pressure_self_cw',
    'pressure_partner_cw',
    'pushing_self_cw',
    'pushing_partner_cw'
  ),
  x_label = c(
    'Received Persuasion',
    'Exerted Persuasion',
    'Received Pressure',
    'Exerted Pressure',
    'Received Plan-Related Pushing',
    'Exerted Plan-Related Pushing'
  ),
  group_var = 'coupleID',
  plot_full_range = TRUE,
  y_limits = c(0, 100),
  y_label = "Same-Day MVPA",
  y_labels = c('Probability of Being Active', 'Minutes of MVPA When Active', 'Overall Expected Minutes of MVPA'),
  , filter_quantiles = .9995
  , font_family = 'Candara'
)
## Warning: Dropping 'draws_df' class as required metadata was removed.
## Warning: Dropping 'draws_df' class as required metadata was removed.
## Warning: Dropping 'draws_df' class as required metadata was removed.
## Warning: Dropping 'draws_df' class as required metadata was removed.
## Warning: Dropping 'draws_df' class as required metadata was removed.
## Warning: Dropping 'draws_df' class as required metadata was removed.
print(conds_eff)

$persuasion_self_cw

## Warning: Removed 88 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning: Removed 43 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Picking joint bandwidth of 0.00821

$persuasion_partner_cw

## Warning: Removed 72 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning: Removed 11 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Picking joint bandwidth of 0.00648

$pressure_self_cw

## Warning: Removed 6 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Picking joint bandwidth of 0.0111

$pressure_partner_cw

## Picking joint bandwidth of 0.0183

$pushing_self_cw

## Picking joint bandwidth of 0.0105

$pushing_partner_cw

## Picking joint bandwidth of 0.0101

Note. This graphic illustrates the relationship between social control and moderate to vigorous physical activity (MVPA) using a Bayesian Hurdle-Lognormal Multilevel Model. The predictor is centered within individuals to examine how deviations from their average social control relate to same-day MVPA. Shaded areas indicate credible intervals, thick lines show fixed effects, and thin lines represent random effects, highlighting variability across couples. The plots display the probability of being active, expected minutes of MVPA when active, and combined predicted MVPA. The bottom density plot visualizes the posterior distributions of slope estimates, transformed to represent multiplicative changes in odds ratios (hurdle component) or expected values. Medians and 95% credible intervals (2.5th and 97.5th percentiles) are shown. Effects are significant, when the 95% credible interval does not overlap 1.

x_label = c(
    'Received Persuasion',
    'Exerted Persuasion',
    'Received Pressure',
    'Exerted Pressure',
    'Received Plan-Related Pushing',
    'Exerted Plan-Related Pushing'
  )

for (i in 1:length(conds_eff)) {
  effname <- names(conds_eff)[i]
  eff_plot <- conds_eff[[i]]
  x_label_i <- x_label[[i]]
  rmarkdown::render(
    "C:/Users/kueng/DataAnalysis/02TimeAndTiesControl/Output/Plots/BeautifulPlotWithNote.Rmd", 
    output_file = paste0('C:/Users/kueng/DataAnalysis/02TimeAndTiesControl/Output/Plots/Graphic_', effname, '.pdf'),
    params = list(
      p_i = eff_plot,
      p_name = effname,
      x_label = x_label_i
      ),
    envir = new.env(),
    quiet = TRUE
  )
}
marginaleffects::avg_predictions(pa_sub)
## 
##  Estimate 2.5 % 97.5 %
##      30.9  29.5   32.4
## 
## Type:  response 
## Columns: estimate, conf.low, conf.high

Device Based MVPA

range(df_double$pa_obj, na.rm = T) 
## [1]   5.75 971.25
hist(df_double$pa_obj, breaks = 50)

df_double$pa_obj_log <- log(df_double$pa_obj)

hist(df_double$pa_obj_log, breaks = 50)

Lognormal Model

formula <- bf(
  pa_obj ~ 
    persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw +
    
    persuasion_self_cb + persuasion_partner_cb +
    pressure_self_cb + pressure_partner_cb +
    pushing_self_cb + pushing_partner_cb +
    
    plan +
    day + weartime_self_cw + weartime_self_cb +
    
    # Random effects
    (persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw | coupleID)
  #, autocor = autocor_str
)



prior1 <- c(
  brms::set_prior("normal(0, 2.5)", class = "b")
  , brms::set_prior("normal(0, 50)", class = "Intercept") 
  , brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
  , brms::set_prior("student_t(3, 0, 2.5)", class = "sigma", lb = 0)
)


#brms::validate_prior(
#  prior1, 
#  formula = formula, 
#  data = df_double, 
#  family = lognormal()
#  )

#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

pa_obj_log <- my_brm(
  mi = use_mi, 
  imputed_data = implist,
  
  formula = formula, 
  prior = prior1,
  data = df_double, 
  family = lognormal(),
  #control = list(adapt_delta = 0.95),
  iter = iterations,
  warmup = warmup,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = file.path("models_cache_brms", paste0("pa_obj_log_gaussian_NOAR", suffix))
)
## Warning: Rows containing NAs were excluded from the model.
pa_obj_log_digest <- digest::digest(pa_obj_log)
check_brms(pa_obj_log, log_pp_check = TRUE)
## 
## Divergences:
## 0 of 40000 iterations ended with a divergence.
## 
## Tree depth:
## 0 of 40000 iterations saturated the maximum tree depth of 10.
## 
## Energy:
## E-BFMI indicated no pathological behavior.

## Using 10 posterior draws for ppc type 'ecdf_overlay' by default.

## Using 10 posterior draws for ppc type 'dens_overlay' by default.

## 
## Computed from 40000 by 3337 log-likelihood matrix.
## 
##          Estimate    SE
## elpd_loo -18777.4  68.9
## p_loo        91.8   4.5
## looic     37554.9 137.8
## ------
## MCSE of elpd_loo is 0.1.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.4, 2.1]).
## 
## All Pareto k estimates are good (k < 0.7).
## See help('pareto-k-diagnostic') for details.
DHARMa.check_brms.all(pa_obj_log, integer = TRUE, outliers_type = 'bootstrap')
## DHARMa:testOutliers with type = binomial may have inflated Type I error rates for integer-valued distributions. To get a more exact result, it is recommended to re-run testOutliers with type = 'bootstrap'. See ?testOutliers for details

## 
##  DHARMa bootstrapped outlier test
## 
## data:  model.check
## outliers at both margin(s) = 24, observations = 3337, p-value < 2.2e-16
## alternative hypothesis: two.sided
##  percent confidence interval:
##  0.001341025 0.005094396
## sample estimates:
## outlier frequency (expected: 0.00314054540005993 ) 
##                                        0.007192089
if (do_priorsense) {
  gc()
  priorsense::powerscale_sensitivity(pa_obj_log, variable = priorsense_vars)
  priorsense::powerscale_plot_dens(pa_obj_log, variable = priorsense_vars)
  priorsense::powerscale_plot_ecdf(pa_obj_log, variable = priorsense_vars)
  priorsense::powerscale_plot_quantities(pa_obj_log, variable = priorsense_vars)
}
# rope range for lognormal model
rope_factor <- sd(log(pa_obj_log$data$pa_obj))
rope_range_log = c(-0.1 * rope_factor, 0.1 * rope_factor)

summarize_brms(
  pa_obj_log, 
  rope_range = rope_range_log,
  model_rows_fixed = model_rows_fixed,
  model_rows_random = model_rows_random,
  model_rownames_fixed = model_rownames_fixed,
  model_rownames_random = model_rownames_random,
  exponentiate = T) %>%
  print_df(rows_to_pack = rows_to_pack)
## Sampling priors, please wait...
## Warning in summarize_brms(pa_obj_log, rope_range = rope_range_log,
## model_rows_fixed = model_rows_fixed, : Coefficients were exponentiated. Double
## check if this was intended.
exp(Est.) SE 95% CI pd ROPE inside ROPE BF BF_Evidence Rhat Bulk_ESS Tail_ESS
Intercept 112.75*** 6.04 [101.44, 125.37] 1.000 [0.94, 1.07] 0.000 >100 Overwhelming Evidence 1.000 5550 12776
Within-Person Effects
Daily persuasion experienced 1.03 0.02 [ 1.00, 1.06] 0.955 [0.94, 1.07] 0.994 0.027 Very Strong Evidence for Null 1.000 23590 29010
Daily persuasion utilized (partner’s view) 1.02 0.02 [ 0.99, 1.05] 0.875 [0.94, 1.07] 0.998 0.012 Very Strong Evidence for Null 1.000 28890 30355
Daily pressure experienced 0.95 0.03 [ 0.88, 1.02] 0.940 [0.94, 1.07] 0.628 0.047 Strong Evidence for Null 1.000 42505 29329
Daily pressure utilized (partner’s view) 0.98 0.03 [ 0.92, 1.05] 0.736 [0.94, 1.07] 0.901 0.016 Very Strong Evidence for Null 1.000 48681 31178
Daily pushing experienced 1.02 0.03 [ 0.96, 1.07] 0.743 [0.94, 1.07] 0.969 0.013 Very Strong Evidence for Null 1.000 31550 27761
Daily pushing utilized (partner’s view) 1.01 0.02 [ 0.97, 1.05] 0.624 [0.94, 1.07] 0.995 0.009 Very Strong Evidence for Null 1.000 43823 28794
Day 0.97 0.03 [ 0.91, 1.04] 0.789 [0.94, 1.07] 0.843 0.019 Very Strong Evidence for Null 1.000 68514 31066
Actionplan 1.09*** 0.03 [ 1.04, 1.14] 1.000 [0.94, 1.07] 0.220 5.363 Moderate Evidence 1.000 57810 30388
Daily weartime 1.00*** 0.00 [ 1.00, 1.00] 1.000 [0.94, 1.07] 1.000 0.091 Strong Evidence for Null 1.000 40853 26910
Between-Person Effects
Mean persuasion experienced 1.10 0.15 [ 0.82, 1.44] 0.750 [0.94, 1.07] 0.280 0.072 Strong Evidence for Null 1.001 5362 9878
Mean persuasion utilized (partner’s view) 0.97 0.14 [ 0.72, 1.28] 0.582 [0.94, 1.07] 0.343 0.057 Strong Evidence for Null 1.001 5445 10506
Mean pressure experienced 1.00 0.14 [ 0.75, 1.34] 0.508 [0.94, 1.07] 0.343 0.058 Strong Evidence for Null 1.001 7720 15395
Mean pressure utilized (partner’s view) 0.97 0.14 [ 0.74, 1.29] 0.579 [0.94, 1.07] 0.345 0.058 Strong Evidence for Null 1.001 7095 14268
Mean pushing experienced 0.95 0.19 [ 0.64, 1.43] 0.598 [0.94, 1.07] 0.240 0.084 Strong Evidence for Null 1.001 8101 14933
Mean pushing utilized (partner’s view) 1.24 0.25 [ 0.83, 1.86] 0.860 [0.94, 1.07] 0.142 0.143 Moderate Evidence for Null 1.001 8132 14838
Mean weartime 1.00 0.00 [ 1.00, 1.00] 0.934 [0.94, 1.07] 1.000 0.000 Very Strong Evidence for Null 1.000 51560 36084
Random Effects
sd(Intercept) 0.30 0.04 [0.23, 0.39] NA NA NA NA NA 1.000 8289 15206
sd(Daily persuasion experienced) 0.05 0.01 [0.02, 0.08] NA NA NA NA NA 1.000 20301 15918
sd(Daily persuasion utilized (partner’s view)) 0.05 0.02 [0.03, 0.09] NA NA NA NA NA 1.000 16365 13439
sd(Daily pressure experienced) 0.04 0.03 [0.00, 0.13] NA NA NA NA NA 1.000 19509 20599
sd(Daily pressure utilized (partner’s view)) 0.03 0.03 [0.00, 0.11] NA NA NA NA NA 1.000 22917 19082
sd(Daily pushing experienced) 0.07 0.04 [0.01, 0.15] NA NA NA NA NA 1.000 9116 8908
sd(Daily pushing utilized (partner’s view)) 0.03 0.03 [0.00, 0.09] NA NA NA NA NA 1.000 15618 19120
Additional Parameters
sigma 0.57 0.01 [0.56, 0.59] NA NA NA NA NA 1.000 66895 28460 
plot(
  bayestestR::p_direction(pa_obj_log),
  priors = TRUE
) + 
  coord_cartesian(xlim = c(-3, 3)) +
  theme_bw()

plot(
  bayestestR::rope(pa_obj_log, range = rope_range_log, ci = 1)
) + theme_bw()
## Possible multicollinearity between b_persuasion_partner_cb and
##   b_persuasion_self_cb (r = 0.89), b_pressure_self_cb and
##   b_persuasion_self_cb (r = 0.73), b_pressure_partner_cb and
##   b_persuasion_self_cb (r = 0.74), b_pressure_self_cb and
##   b_persuasion_partner_cb (r = 0.71), b_pressure_partner_cb and
##   b_persuasion_partner_cb (r = 0.77), b_pushing_partner_cb and
##   b_pushing_self_cb (r = 0.82). This might lead to inappropriate results.
##   See 'Details' in '?rope'.

# Nothing significant, no plots

Affect

range(df_double$aff, na.rm = T)
## [1] 0 5
hist(df_double$aff, breaks = 15)

Gaussian

formula <- bf(
  aff ~ 
    persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw +
    
    persuasion_self_cb + persuasion_partner_cb +
    pressure_self_cb + pressure_partner_cb +
    pushing_self_cb + pushing_partner_cb +
    
    plan +
    day +
    
    # Random effects
    (persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw | coupleID)
  #, autocor = autocor_str
)


prior1 <- c(
  brms::set_prior("normal(0, 5)", class = "b")
  ,brms::set_prior("normal(0, 20)", class = "Intercept", lb=1, ub=6)
  , brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
  , brms::set_prior("student_t(3, 0, 2.5)", class = "sigma", lb = 0)
)

#brms::validate_prior(
#  prior1, 
#  formula = formula, 
#  data = df_double, 
#  family = gaussian()
#  )

#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

mood_gauss <- my_brm(
  mi = use_mi, 
  imputed_data = implist,
  
  formula = formula, 
  prior = prior1,
  data = df_double, 
  family = gaussian(),
  #control = list(adapt_delta = 0.95, max_treedepth = 15),
  iter = iterations,
  warmup = warmup,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = file.path("models_cache_brms", paste0("mood_gauss_NOAR", suffix))
)
## Warning: Rows containing NAs were excluded from the model.
mood_gauss_digest <- digest::digest(mood_gauss)
check_brms(mood_gauss, log_pp_check = FALSE)
## 
## Divergences:
## 0 of 40000 iterations ended with a divergence.
## 
## Tree depth:
## 0 of 40000 iterations saturated the maximum tree depth of 10.
## 
## Energy:
## E-BFMI indicated no pathological behavior.

## Using 10 posterior draws for ppc type 'ecdf_overlay' by default.

## Using 10 posterior draws for ppc type 'dens_overlay' by default.

## Warning: Found 1 observations with a pareto_k > 0.7 in model 'model'. We
## recommend to set 'moment_match = TRUE' in order to perform moment matching for
## problematic observations.
## 
## Computed from 40000 by 3736 log-likelihood matrix.
## 
##          Estimate    SE
## elpd_loo  -5184.4  59.1
## p_loo        75.9   3.3
## looic     10368.9 118.2
## ------
## MCSE of elpd_loo is NA.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.5, 2.1]).
## 
## Pareto k diagnostic values:
##                          Count Pct.    Min. ESS
## (-Inf, 0.7]   (good)     3735  100.0%  5669    
##    (0.7, 1]   (bad)         1    0.0%  <NA>    
##    (1, Inf)   (very bad)    0    0.0%  <NA>    
## See help('pareto-k-diagnostic') for details.
DHARMa.check_brms.all(mood_gauss, integer = FALSE)

## 
##  DHARMa outlier test based on exact binomial test with approximate
##  expectations
## 
## data:  model.check
## outliers at both margin(s) = 27, observations = 3736, p-value =
## 2.533e-08
## alternative hypothesis: true probability of success is not equal to 0.001998002
## 95 percent confidence interval:
##  0.004767871 0.010497581
## sample estimates:
## frequency of outliers (expected: 0.001998001998002 ) 
##                                          0.007226981
if (do_priorsense) {
  gc()
  priorsense::powerscale_sensitivity(mood_gauss, variable = priorsense_vars)
  priorsense::powerscale_plot_dens(mood_gauss, variable = priorsense_vars)
  priorsense::powerscale_plot_ecdf(mood_gauss, variable = priorsense_vars)
  priorsense::powerscale_plot_quantities(mood_gauss, variable = priorsense_vars)
}
summarize_brms(
  mood_gauss, 
  model_rows_fixed = model_rows_fixed,
  model_rows_random = model_rows_random,
  model_rownames_fixed = model_rownames_fixed,
  model_rownames_random = model_rownames_random,
  exponentiate = F) %>%
  print_df(rows_to_pack = rows_to_pack)
## Sampling priors, please wait...
Est. SE 95% CI pd ROPE inside ROPE BF BF_Evidence Rhat Bulk_ESS Tail_ESS
Intercept 3.66*** 0.10 [ 3.45, 3.87] 1.000 [-0.11, 0.11] 0.000 >100 Overwhelming Evidence 1.001 4553 9173
Within-Person Effects
Daily persuasion experienced 0.00 0.02 [-0.04, 0.04] 0.519 [-0.11, 0.11] 1.000 0.004 Very Strong Evidence for Null 1.000 37309 29134
Daily persuasion utilized (partner’s view) 0.02 0.02 [-0.03, 0.07] 0.808 [-0.11, 0.11] 1.000 0.007 Very Strong Evidence for Null 1.000 30062 28707
Daily pressure experienced -0.03 0.05 [-0.14, 0.07] 0.724 [-0.11, 0.11] 0.931 0.012 Very Strong Evidence for Null 1.000 40161 28685
Daily pressure utilized (partner’s view) -0.02 0.05 [-0.14, 0.08] 0.675 [-0.11, 0.11] 0.933 0.012 Very Strong Evidence for Null 1.000 33036 22902
Daily pushing experienced 0.00 0.03 [-0.06, 0.07] 0.561 [-0.11, 0.11] 0.999 0.007 Very Strong Evidence for Null 1.000 39938 30978
Daily pushing utilized (partner’s view) 0.07 0.03 [ 0.00, 0.13] 0.972 [-0.11, 0.11] 0.917 0.047 Strong Evidence for Null 1.000 32350 26848
Day 0.26*** 0.06 [ 0.15, 0.36] 1.000 [-0.11, 0.11] 0.004 >100 Overwhelming Evidence 1.000 63577 31006
Actionplan 0.09* 0.04 [ 0.02, 0.16] 0.994 [-0.11, 0.11] 0.753 0.170 Moderate Evidence for Null 1.000 53515 29497
Daily weartime NA NA NA NA NA NA NA NA NA NA NA
Between-Person Effects
Mean persuasion experienced 0.34 0.28 [-0.21, 0.89] 0.889 [-0.11, 0.11] 0.157 0.116 Moderate Evidence for Null 1.001 4317 8353
Mean persuasion utilized (partner’s view) 0.23 0.28 [-0.33, 0.77] 0.792 [-0.11, 0.11] 0.234 0.077 Strong Evidence for Null 1.001 4313 8549
Mean pressure experienced -0.30 0.27 [-0.84, 0.24] 0.864 [-0.11, 0.11] 0.186 0.100 Strong Evidence for Null 1.001 5213 10323
Mean pressure utilized (partner’s view) -0.31 0.27 [-0.85, 0.23] 0.871 [-0.11, 0.11] 0.179 0.102 Moderate Evidence for Null 1.001 5234 9935
Mean pushing experienced 0.19 0.38 [-0.57, 0.95] 0.692 [-0.11, 0.11] 0.210 0.085 Strong Evidence for Null 1.001 6886 13871
Mean pushing utilized (partner’s view) 0.34 0.38 [-0.42, 1.10] 0.820 [-0.11, 0.11] 0.160 0.116 Moderate Evidence for Null 1.001 6875 13862
Mean weartime NA NA NA NA NA NA NA NA NA NA NA
Random Effects
sd(Intercept) 0.60 0.07 [0.47, 0.78] NA NA NA NA NA 1.001 7034 14642
sd(Daily persuasion experienced) 0.04 0.03 [0.00, 0.10] NA NA NA NA NA 1.001 10094 13693
sd(Daily persuasion utilized (partner’s view)) 0.07 0.03 [0.01, 0.13] NA NA NA NA NA 1.000 10348 8520
sd(Daily pressure experienced) 0.07 0.06 [0.00, 0.24] NA NA NA NA NA 1.000 16978 18579
sd(Daily pressure utilized (partner’s view)) 0.08 0.07 [0.00, 0.26] NA NA NA NA NA 1.000 15464 18906
sd(Daily pushing experienced) 0.05 0.04 [0.00, 0.14] NA NA NA NA NA 1.000 15219 18471
sd(Daily pushing utilized (partner’s view)) 0.07 0.05 [0.00, 0.17] NA NA NA NA NA 1.000 13502 14099
Additional Parameters
sigma 0.96 0.01 [0.94, 0.98] NA NA NA NA NA 1.000 58119 30453 
plot(
  bayestestR::p_direction(mood_gauss),
  priors = TRUE
)  + 
  coord_cartesian(xlim = c(-3, 3)) +
  theme_bw()
## Warning in `==.default`(dens$Parameter, parameter): longer object length is not
## a multiple of shorter object length
## Warning in is.na(e1) | is.na(e2): longer object length is not a multiple of
## shorter object length

plot(
  bayestestR::rope(mood_gauss, ci = 1)
) + theme_bw()
## Possible multicollinearity between b_pressure_self_cb and
##   b_persuasion_self_cb (r = 0.83), b_pressure_partner_cb and
##   b_persuasion_self_cb (r = 0.8), b_pressure_self_cb and
##   b_persuasion_partner_cb (r = 0.8), b_pressure_partner_cb and
##   b_persuasion_partner_cb (r = 0.82), b_pressure_partner_cb and
##   b_pressure_self_cb (r = 0.76), b_pushing_partner_cb and
##   b_pushing_self_cb (r = 0.88). This might lead to inappropriate results.
##   See 'Details' in '?rope'.

conditional_spaghetti(
  mood_gauss, 
  effects = c('pushing_partner_cw'),
  group_var = 'coupleID',
  plot_full_range = TRUE
)

$pushing_partner_cw 

marginaleffects::avg_predictions(mood_gauss)
## 
##  Estimate 2.5 % 97.5 %
##      3.82  3.79   3.85
## 
## Type:  response 
## Columns: estimate, conf.low, conf.high

Reactance

range(df_double$reactance, na.rm = T) 
## [1] 0 5
hist(df_double$reactance, breaks = 7) 

hist(log(df_double$reactance+0.1), breaks = 10)

Ordinal

df_double$reactance_ordinal <- factor(df_double$reactance,
                                      levels = 0:5, 
                                      ordered = TRUE)

formula <- bf(
  reactance_ordinal ~ 
    persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw +
    
    persuasion_self_cb + persuasion_partner_cb +
    pressure_self_cb + pressure_partner_cb +
    pushing_self_cb + pushing_partner_cb +
    
    plan +
    day +
    
    # Random effects
    (persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw | coupleID)
  #, autocor = autocor_str
)


prior1 <- c(
  brms::set_prior("normal(0, 2.5)", class = "b")
  , brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
)


#brms::validate_prior(
#  prior1, 
#  formula = formula, 
#  data = df_double, 
#  family = cumulative() # HURDLE_CUMULATIVE
#  )


#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

reactance_ordinal <- my_brm(
  mi = use_mi, 
  imputed_data = implist,
  
  formula = formula, 
  prior = prior1,
  data = df_double, 
  family = brms::cumulative(),
  #control = list(adapt_delta = 0.95),
  iter = iterations,
  warmup = warmup,
  chains = 4,
  cores = 4,
  seed = 7777
  , file = file.path("models_cache_brms", paste0("reactance_ordinal_NOARNOAR_", suffix))
)
## Warning: Rows containing NAs were excluded from the model.
reactance_ordinal_digest <- digest::digest(reactance_ordinal)
check_brms(reactance_ordinal)
## 
## Divergences:
## 0 of 40000 iterations ended with a divergence.
## 
## Tree depth:
## 0 of 40000 iterations saturated the maximum tree depth of 10.
## 
## Energy:
## E-BFMI indicated no pathological behavior.

## Using 10 posterior draws for ppc type 'ecdf_overlay' by default.

## Using 10 posterior draws for ppc type 'dens_overlay' by default.

## Warning: Found 6 observations with a pareto_k > 0.7 in model 'model'. We
## recommend to set 'moment_match = TRUE' in order to perform moment matching for
## problematic observations.
## 
## Computed from 40000 by 756 log-likelihood matrix.
## 
##          Estimate   SE
## elpd_loo   -682.8 32.0
## p_loo        75.0  5.5
## looic      1365.6 64.0
## ------
## MCSE of elpd_loo is NA.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.4, 1.7]).
## 
## Pareto k diagnostic values:
##                          Count Pct.    Min. ESS
## (-Inf, 0.7]   (good)     750   99.2%   431     
##    (0.7, 1]   (bad)        6    0.8%   <NA>    
##    (1, Inf)   (very bad)   0    0.0%   <NA>    
## See help('pareto-k-diagnostic') for details.
DHARMa.check_brms.all(reactance_ordinal, outliers_type = 'bootstrap')

## 
##  DHARMa bootstrapped outlier test
## 
## data:  model.check
## outliers at both margin(s) = 2, observations = 756, p-value = 0.08
## alternative hypothesis: two.sided
##  percent confidence interval:
##  0.000000000 0.002645503
## sample estimates:
## outlier frequency (expected: 0.000343915343915344 ) 
##                                         0.002645503
if (do_priorsense) {
  gc()
  priorsense::powerscale_sensitivity(reactance_ordinal, variable = priorsense_vars)
  priorsense::powerscale_plot_dens(reactance_ordinal, variable = priorsense_vars)
  priorsense::powerscale_plot_ecdf(reactance_ordinal, variable = priorsense_vars)
  priorsense::powerscale_plot_quantities(reactance_ordinal, variable = priorsense_vars)
}
summarize_brms(
  reactance_ordinal, 
  rope_range = c(-0.18, 0.18),
  model_rows_fixed = model_rows_fixed_ordinal,
  model_rows_random = model_rows_random_ordinal,
  model_rownames_fixed = model_rownames_fixed_ordinal,
  model_rownames_random = model_rownames_random_ordinal,
  exponentiate = T) %>%
  print_df(rows_to_pack = rows_to_pack_ordinal)
## Sampling priors, please wait...
OR SE 95% CI pd ROPE inside ROPE BF BF_Evidence Rhat Bulk_ESS Tail_ESS
Intercepts
Intercept NA NA NA NA NA NA NA NA NA NA NA
Intercept[1] 3.44*** 1.02 [ 1.94, 6.27] 1.000 [0.84, 1.20] 0.000 >100 Overwhelming Evidence 1.000 34475 33538
Intercept[2] 7.48*** 2.32 [ 4.13, 13.93] 1.000 [0.84, 1.20] 0.000 >100 Overwhelming Evidence 1.000 34550 32342
Intercept[3] 20.85*** 6.89 [ 11.02, 41.02] 1.000 [0.84, 1.20] 0.000 >100 Overwhelming Evidence 1.000 36864 31642
Intercept[4] 91.11*** 34.90 [ 44.35, 198.17] 1.000 [0.84, 1.20] 0.000 >100 Overwhelming Evidence 1.000 41991 33122
Intercept[5] 3062.08*** 2032.71 [949.61, 12641.71] 1.000 [0.84, 1.20] 0.000 >100 Overwhelming Evidence 1.000 53681 30726
Within-Person Effects
Daily persuasion experienced 0.84* 0.07 [ 0.71, 0.99] 0.981 [0.84, 1.20] 0.553 0.288 Moderate Evidence for Null 1.000 39334 29508
Daily persuasion utilized (partner’s view) 1.02 0.10 [ 0.83, 1.23] 0.590 [0.84, 1.20] 0.925 0.040 Strong Evidence for Null 1.000 36576 30022
Daily pressure experienced 1.84* 0.36 [ 1.16, 2.65] 0.993 [0.84, 1.20] 0.031 2.516 Weak Evidence 1.000 21421 23286
Daily pressure utilized (partner’s view) 1.24 0.29 [ 0.70, 2.14] 0.817 [0.84, 1.20] 0.366 0.162 Moderate Evidence for Null 1.000 26173 22280
Daily pushing experienced 1.21 0.13 [ 0.98, 1.51] 0.962 [0.84, 1.20] 0.449 0.220 Moderate Evidence for Null 1.000 31297 29838
Daily pushing utilized (partner’s view) 0.92 0.11 [ 0.71, 1.19] 0.741 [0.84, 1.20] 0.769 0.062 Strong Evidence for Null 1.000 37378 28440
Day 1.49 0.51 [ 0.76, 2.90] 0.876 [0.84, 1.20] 0.219 0.269 Moderate Evidence for Null 1.000 50318 30070
Actionplan 0.82 0.21 [ 0.50, 1.34] 0.788 [0.84, 1.20] 0.399 0.137 Moderate Evidence for Null 1.000 47755 33168
Daily weartime NA NA NA NA NA NA NA NA NA NA NA
Between-Person Effects
Mean persuasion experienced 1.09 0.56 [ 0.39, 3.05] 0.566 [0.84, 1.20] 0.272 0.203 Moderate Evidence for Null 1.001 16726 24563
Mean persuasion utilized (partner’s view) 1.39 0.79 [ 0.46, 4.37] 0.719 [0.84, 1.20] 0.214 0.262 Moderate Evidence for Null 1.000 17310 24106
Mean pressure experienced 3.51* 1.92 [ 1.20, 10.46] 0.989 [0.84, 1.20] 0.020 3.094 Moderate Evidence 1.001 18478 25633
Mean pressure utilized (partner’s view) 1.16 0.65 [ 0.36, 3.48] 0.602 [0.84, 1.20] 0.242 0.232 Moderate Evidence for Null 1.000 18525 24890
Mean pushing experienced 1.26 0.94 [ 0.29, 5.83] 0.621 [0.84, 1.20] 0.183 0.313 Weak Evidence for Null 1.000 22187 26801
Mean pushing utilized (partner’s view) 0.11* 0.10 [ 0.02, 0.63] 0.994 [0.84, 1.20] 0.008 8.118 Moderate Evidence 1.000 27918 29396
Mean weartime NA NA NA NA NA NA NA NA NA NA NA
Random Effects
sd(Intercept) 0.80 0.20 [0.47, 1.27] NA NA NA NA NA 1.000 17242 26661
sd(Daily persuasion experienced) 0.17 0.12 [0.01, 0.43] NA NA NA NA NA 1.001 7879 13881
sd(Daily persuasion utilized (partner’s view)) 0.21 0.14 [0.01, 0.52] NA NA NA NA NA 1.000 11746 14888
sd(Daily pressure experienced) 0.55 0.25 [0.09, 1.16] NA NA NA NA NA 1.001 9564 8028
sd(Daily pressure utilized (partner’s view)) 0.42 0.39 [0.02, 1.58] NA NA NA NA NA 1.000 10080 18938
sd(Daily pushing experienced) 0.20 0.13 [0.01, 0.50] NA NA NA NA NA 1.000 12311 15010
sd(Daily pushing utilized (partner’s view)) 0.15 0.14 [0.01, 0.61] NA NA NA NA NA 1.000 17143 21955
Additional Parameters
sigma NA NA NA NA NA NA NA NA NA NA NA
disc 1.00 0.00 [1.00, 1.00] NA NA NA NA NA NA NA NA 
plot(
  bayestestR::p_direction(reactance_ordinal),
  priors = TRUE
) + 
  coord_cartesian(xlim = c(-6, 6)) +
  theme_bw()
## Warning in `==.default`(dens$Parameter, parameter): longer object length is not
## a multiple of shorter object length
## Warning in is.na(e1) | is.na(e2): longer object length is not a multiple of
## shorter object length

plot(
  bayestestR::rope(reactance_ordinal, range = c(-0.18, 0.18), ci = 1)
) + theme_bw()
## Possible multicollinearity between b_Intercept[4] and b_Intercept[2] (r
##   = 0.8), b_Intercept[4] and b_Intercept[3] (r = 0.86), b_pressure_self_cb
##   and b_persuasion_self_cb (r = 0.72), b_pressure_partner_cb and
##   b_persuasion_partner_cb (r = 0.79). This might lead to inappropriate
##   results. See 'Details' in '?rope'.

conditional_spaghetti(
  reactance_ordinal, 
  effects = c('persuasion_self_cw', 'pressure_self_cw')
  , group_var = 'coupleID'
  #, n_groups = 15
  , plot_full_range = T
)

\(persuasion_self_cw <img src="01_FinalModels_files/figure-html/report_reactance_ordinal-3.png" width="2400" />\)pressure_self_cw 

marginaleffects::avg_predictions(reactance_ordinal)
## 
##  Group Estimate   2.5 % 97.5 %
##      0  0.68198 0.65457 0.7080
##      1  0.09409 0.07478 0.1155
##      2  0.08360 0.06612 0.1039
##      3  0.06980 0.05432 0.0878
##      4  0.06387 0.05134 0.0776
##      5  0.00527 0.00171 0.0116
## 
## Type:  response 
## Columns: group, estimate, conf.low, conf.high

Binary

introduce_binary_reactance <- function(data) {
  data$is_reactance <- factor(data$reactance > 0, levels = c(FALSE, TRUE), labels = c(0, 1))
  return(data)
}



df_double <- introduce_binary_reactance(df_double)
if (use_mi) {
  for (i in seq_along(implist)) {
    implist[[i]] <- introduce_binary_reactance(implist[[i]])
  }
}


formula <- bf(
  is_reactance ~ 
    persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw +
    
    persuasion_self_cb + persuasion_partner_cb +
    pressure_self_cb + pressure_partner_cb +
    pushing_self_cb + pushing_partner_cb +
    
    plan +
    day +
    
    # Random effects
    (persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw | coupleID)
  #, autocor = autocor_str
  )



prior1 <- c(
  brms::set_prior("normal(0, 2.5)", class = "b")
  , brms::set_prior("normal(0, 10)", class = "Intercept", lb=0, ub=5) 
  , brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
)


#brms::validate_prior(
#  prior1, 
#  formula = formula, 
#  data = df_double, 
#  family = bernoulli()
#  )



#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

is_reactance <- my_brm(
  mi = use_mi, 
  imputed_data = implist,
  
  formula = formula, 
  prior = prior1,
  data = df_double, 
  family = brms::bernoulli(),
  #control = list(adapt_delta = 0.95, max_treedepth = 15),
  iter = iterations,
  warmup = warmup,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = file.path("models_cache_brms", paste0("is_reactance_NOARNOAR_", suffix))
  #, file_refit = 'always'
)
## Warning: Rows containing NAs were excluded from the model.
is_reactance_digest <- digest::digest(is_reactance)
check_brms(is_reactance)
## 
## Divergences:
## 0 of 40000 iterations ended with a divergence.
## 
## Tree depth:
## 0 of 40000 iterations saturated the maximum tree depth of 10.
## 
## Energy:
## E-BFMI indicated no pathological behavior.

## Using 10 posterior draws for ppc type 'ecdf_overlay' by default.

## Using 10 posterior draws for ppc type 'dens_overlay' by default.

## Warning: Found 41 observations with a pareto_k > 0.7 in model 'model'. We
## recommend to set 'moment_match = TRUE' in order to perform moment matching for
## problematic observations.
## 
## Computed from 40000 by 756 log-likelihood matrix.
## 
##          Estimate   SE
## elpd_loo   -364.3 16.0
## p_loo        81.8  6.0
## looic       728.6 32.0
## ------
## MCSE of elpd_loo is NA.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.3, 1.4]).
## 
## Pareto k diagnostic values:
##                          Count Pct.    Min. ESS
## (-Inf, 0.7]   (good)     715   94.6%   684     
##    (0.7, 1]   (bad)       40    5.3%   <NA>    
##    (1, Inf)   (very bad)   1    0.1%   <NA>    
## See help('pareto-k-diagnostic') for details.
DHARMa.check_brms.all(is_reactance, integer = FALSE)

## 
##  DHARMa outlier test based on exact binomial test with approximate
##  expectations
## 
## data:  model.check
## outliers at both margin(s) = 1, observations = 756, p-value = 1
## alternative hypothesis: true probability of success is not equal to 0.001998002
## 95 percent confidence interval:
##  0.0000334886 0.0073476538
## sample estimates:
## frequency of outliers (expected: 0.001998001998002 ) 
##                                          0.001322751
if (do_priorsense) {
  gc()
  priorsense::powerscale_sensitivity(is_reactance, variable = priorsense_vars)
  priorsense::powerscale_plot_dens(is_reactance, variable = priorsense_vars)
  priorsense::powerscale_plot_ecdf(is_reactance, variable = priorsense_vars)
  priorsense::powerscale_plot_quantities(is_reactance, variable = priorsense_vars)
}
summarize_brms(
  is_reactance, 
  model_rows_fixed = model_rows_fixed,
  model_rows_random = model_rows_random,
  model_rownames_fixed = model_rownames_fixed,
  model_rownames_random = model_rownames_random,
  exponentiate = T) %>%
  print_df(rows_to_pack = rows_to_pack)
## Sampling priors, please wait...
OR SE 95% CI pd ROPE inside ROPE BF BF_Evidence Rhat Bulk_ESS Tail_ESS
Intercept 0.33*** 0.11 [0.17, 0.61] 1.000 [0.83, 1.20] 0.002 22.786 Strong Evidence 1.000 27240 29839
Within-Person Effects
Daily persuasion experienced 0.84 0.08 [0.68, 1.02] 0.965 [0.83, 1.20] 0.533 0.206 Moderate Evidence for Null 1.000 23994 22476
Daily persuasion utilized (partner’s view) 1.11 0.16 [0.84, 1.52] 0.777 [0.83, 1.20] 0.666 0.078 Strong Evidence for Null 1.000 19267 22853
Daily pressure experienced 1.97* 0.64 [1.01, 4.51] 0.977 [0.83, 1.20] 0.054 1.283 Weak Evidence 1.000 16689 17375
Daily pressure utilized (partner’s view) 1.44 0.61 [0.58, 4.26] 0.811 [0.83, 1.20] 0.227 0.249 Moderate Evidence for Null 1.000 17194 17009
Daily pushing experienced 1.32* 0.17 [1.03, 1.75] 0.985 [0.83, 1.20] 0.220 0.572 Weak Evidence for Null 1.000 23616 24332
Daily pushing utilized (partner’s view) 0.90 0.17 [0.61, 1.34] 0.702 [0.83, 1.20] 0.588 0.087 Strong Evidence for Null 1.000 24981 24147
Day 1.66 0.65 [0.77, 3.61] 0.902 [0.83, 1.20] 0.165 0.360 Weak Evidence for Null 1.000 36962 29582
Actionplan 0.81 0.23 [0.47, 1.40] 0.779 [0.83, 1.20] 0.375 0.146 Moderate Evidence for Null 1.000 34447 29971
Daily weartime NA NA NA NA NA NA NA NA NA NA NA
Between-Person Effects
Mean persuasion experienced 1.89 1.14 [0.58, 6.55] 0.854 [0.83, 1.20] 0.137 0.411 Weak Evidence for Null 1.000 12950 21238
Mean persuasion utilized (partner’s view) 1.84 1.21 [0.50, 7.09] 0.825 [0.83, 1.20] 0.141 0.402 Weak Evidence for Null 1.000 14420 22244
Mean pressure experienced 18.13** 19.21 [2.50, 165.12] 0.999 [0.83, 1.20] 0.002 31.950 Very Strong Evidence 1.000 15560 21627
Mean pressure utilized (partner’s view) 2.34 2.57 [0.25, 19.20] 0.778 [0.83, 1.20] 0.099 0.599 Weak Evidence for Null 1.000 16606 23407
Mean pushing experienced 0.87 0.87 [0.12, 6.33] 0.557 [0.83, 1.20] 0.143 0.397 Weak Evidence for Null 1.000 16461 24158
Mean pushing utilized (partner’s view) 0.08* 0.09 [0.01, 0.65] 0.991 [0.83, 1.20] 0.008 7.125 Moderate Evidence 1.000 19649 26087
Mean weartime NA NA NA NA NA NA NA NA NA NA NA
Random Effects
sd(Intercept) 1.16 0.25 [0.74, 1.74] NA NA NA NA NA 1.000 12403 21667
sd(Daily persuasion experienced) 0.21 0.14 [0.01, 0.52] NA NA NA NA NA 1.001 6876 11478
sd(Daily persuasion utilized (partner’s view)) 0.48 0.20 [0.12, 0.98] NA NA NA NA NA 1.001 10749 9966
sd(Daily pressure experienced) 1.08 0.56 [0.13, 2.46] NA NA NA NA NA 1.000 7167 6747
sd(Daily pressure utilized (partner’s view)) 0.83 0.67 [0.04, 2.75] NA NA NA NA NA 1.000 11499 14462
sd(Daily pushing experienced) 0.25 0.16 [0.02, 0.61] NA NA NA NA NA 1.000 11923 13576
sd(Daily pushing utilized (partner’s view)) 0.26 0.23 [0.01, 0.96] NA NA NA NA NA 1.000 12936 17977
Additional Parameters
sigma NA NA NA NA NA NA NA NA NA NA NA 
plot(
  bayestestR::p_direction(is_reactance),
  priors = TRUE
) + 
  coord_cartesian(xlim = c(-6, 6)) +
  theme_bw()
## Warning in `==.default`(dens$Parameter, parameter): longer object length is not
## a multiple of shorter object length
## Warning in is.na(e1) | is.na(e2): longer object length is not a multiple of
## shorter object length

plot(
  bayestestR::rope(is_reactance, ci = 1)
) + theme_bw()

conditional_spaghetti(
  is_reactance, 
  effects = c('pressure_self_cw', 'pushing_self_cw'),
  group_var = 'coupleID',
  plot_full_range = TRUE
)

\(pressure_self_cw <img src="01_FinalModels_files/figure-html/report_is_reactance-3.png" width="2400" />\)pushing_self_cw 

marginaleffects::avg_predictions(is_reactance)
## 
##  Estimate 2.5 % 97.5 %
##     0.329 0.304  0.354
## 
## Type:  response 
## Columns: estimate, conf.low, conf.high
hypothesis(is_reactance, "pressure_self_cw > pushing_self_cw")
## Hypothesis Tests for class b:
##                 Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio
## 1 (pressure_self_cw... > 0     0.42      0.41    -0.22     1.11       6.45
##   Post.Prob Star
## 1      0.87     
## ---
## 'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
## '*': For one-sided hypotheses, the posterior probability exceeds 95%;
## for two-sided hypotheses, the value tested against lies outside the 95%-CI.
## Posterior probabilities of point hypotheses assume equal prior probabilities.

Report All Models

process_model_component <- function(obj) {
  # Convert to string, modify, and evaluate
  name <- deparse(substitute(obj))
  if (report_hurdle) name <- paste0(name, '_hu')
  if (report_ordinal) name <- paste0(name, '_ordinal')
  return(get(name, envir = parent.frame()))
}



all_models <- report_side_by_side(
  pa_sub,
  pa_obj_log,
  mood_gauss,
  reactance_ordinal,
  is_reactance,
  
  stats_to_report = c('CI', 'pd'),
  
  model_rows_random = process_model_component(model_rows_random),
  model_rows_fixed = process_model_component(model_rows_fixed),
  model_rownames_random = process_model_component(model_rownames_random),
  model_rownames_fixed = process_model_component(model_rownames_fixed)
)

[1] “pa_sub”

## Warning in summarize_brms(model, exponentiate = exponentiate, stats_to_report =
## stats_to_report, : Coefficients were exponentiated. Double check if this was
## intended.

[1] “pa_obj_log”

## Warning in summarize_brms(model, exponentiate = exponentiate, stats_to_report =
## stats_to_report, : Coefficients were exponentiated. Double check if this was
## intended.

[1] “mood_gauss” [1] “reactance_ordinal” [1] “is_reactance”

# pretty printing

summary_all_models <- all_models %>%
  print_df(rows_to_pack = process_model_component(rows_to_pack)) %>%
  add_header_above(
    c(" ", "Subjective MVPA Hurdle Lognormal" = (length(all_models) / 5),  
      "Device-Based MVPA Log (Gaussian)" = (length(all_models) / 5), 
      "Mood Gaussian" = (length(all_models) / 5),
      "Reactance Ordinal" = (length(all_models) / 5),
      "Reactance Dichotome" = (length(all_models) / 5))
  )


export_xlsx(
  summary_all_models, 
  rows_to_pack = process_model_component(rows_to_pack),
  file.path("Output", "AllModelsFinal.xlsx"), 
  merge_option = 'header', 
  simplify_2nd_row = TRUE,
  line_above_rows = c(1,2),
  line_below_rows = c(-1)
)

  
print(summary_all_models)
Subjective MVPA Hurdle Lognormal
Device-Based MVPA Log (Gaussian)
Mood Gaussian
Reactance Ordinal
Reactance Dichotome
exp(Est.) pa_sub 95% CI pa_sub pd pa_sub exp(Est.) pa_obj_log 95% CI pa_obj_log pd pa_obj_log Est. mood_gauss 95% CI mood_gauss pd mood_gauss OR reactance_ordinal 95% CI reactance_ordinal pd reactance_ordinal OR is_reactance 95% CI is_reactance pd is_reactance
Intercept 38.88*** [33.95, 44.46] 1.000 112.75*** [101.44, 125.37] 1.000 3.66*** [ 3.45, 3.87] 1.000 NA NA NA 0.33*** [0.17, 0.61] 1.000
Hurdle Intercept 0.27*** [ 0.20, 0.36] 1.000 NA NA NA NA NA NA NA NA NA NA NA NA
Conditional Within-Person Effects
Daily persuasion experienced 1.03 [ 0.98, 1.09] 0.862 1.03 [ 1.00, 1.06] 0.955 0.00 [-0.04, 0.04] 0.519 0.84* [ 0.71, 0.99] 0.981 0.84 [0.68, 1.02] 0.965
Daily persuasion utilized (partner’s view) 1.03 [ 0.99, 1.08] 0.922 1.02 [ 0.99, 1.05] 0.875 0.02 [-0.03, 0.07] 0.808 1.02 [ 0.83, 1.23] 0.590 1.11 [0.84, 1.52] 0.777
Daily pressure experienced 0.91 [ 0.82, 1.00] 0.972 0.95 [ 0.88, 1.02] 0.940 -0.03 [-0.14, 0.07] 0.724 1.84* [ 1.16, 2.65] 0.993 1.97* [1.01, 4.51] 0.977
Daily pressure utilized (partner’s view) 0.94 [ 0.86, 1.03] 0.917 0.98 [ 0.92, 1.05] 0.736 -0.02 [-0.14, 0.08] 0.675 1.24 [ 0.70, 2.14] 0.817 1.44 [0.58, 4.26] 0.811
Daily pushing experienced 0.99 [ 0.93, 1.06] 0.614 1.02 [ 0.96, 1.07] 0.743 0.00 [-0.06, 0.07] 0.561 1.21 [ 0.98, 1.51] 0.962 1.32* [1.03, 1.75] 0.985
Daily pushing utilized (partner’s view) 0.97 [ 0.91, 1.03] 0.859 1.01 [ 0.97, 1.05] 0.624 0.07 [ 0.00, 0.13] 0.972 0.92 [ 0.71, 1.19] 0.741 0.90 [0.61, 1.34] 0.702
Day 0.99 [ 0.88, 1.12] 0.559 0.97 [ 0.91, 1.04] 0.789 0.26*** [ 0.15, 0.36] 1.000 1.49 [ 0.76, 2.90] 0.876 1.66 [0.77, 3.61] 0.902
Actionplan 1.36*** [ 1.25, 1.49] 1.000 1.09*** [ 1.04, 1.14] 1.000 0.09* [ 0.02, 0.16] 0.994 0.82 [ 0.50, 1.34] 0.788 0.81 [0.47, 1.40] 0.779
Daily weartime NA NA NA 1.00*** [ 1.00, 1.00] 1.000 NA NA NA NA NA NA NA NA NA
Conditional Between-Person Effects
Mean persuasion experienced 1.04 [ 0.76, 1.40] 0.588 1.10 [ 0.82, 1.44] 0.750 0.34 [-0.21, 0.89] 0.889 1.09 [ 0.39, 3.05] 0.566 1.89 [0.58, 6.55] 0.854
Mean persuasion utilized (partner’s view) 0.98 [ 0.73, 1.33] 0.543 0.97 [ 0.72, 1.28] 0.582 0.23 [-0.33, 0.77] 0.792 1.39 [ 0.46, 4.37] 0.719 1.84 [0.50, 7.09] 0.825
Mean pressure experienced 1.18 [ 0.84, 1.69] 0.831 1.00 [ 0.75, 1.34] 0.508 -0.30 [-0.84, 0.24] 0.864 3.51* [ 1.20, 10.46] 0.989 18.13** [2.50, 165.12] 0.999
Mean pressure utilized (partner’s view) 0.92 [ 0.64, 1.32] 0.678 0.97 [ 0.74, 1.29] 0.579 -0.31 [-0.85, 0.23] 0.871 1.16 [ 0.36, 3.48] 0.602 2.34 [0.25, 19.20] 0.778
Mean pushing experienced 1.21 [ 0.78, 1.91] 0.800 0.95 [ 0.64, 1.43] 0.598 0.19 [-0.57, 0.95] 0.692 1.26 [ 0.29, 5.83] 0.621 0.87 [0.12, 6.33] 0.557
Mean pushing utilized (partner’s view) 1.33 [ 0.84, 2.11] 0.893 1.24 [ 0.83, 1.86] 0.860 0.34 [-0.42, 1.10] 0.820 0.11* [ 0.02, 0.63] 0.994 0.08* [0.01, 0.65] 0.991
Mean weartime NA NA NA 1.00 [ 1.00, 1.00] 0.934 NA NA NA NA NA NA NA NA NA
Hurdle Within-Person Effects
Hu Daily persuasion experienced 1.59*** [ 1.39, 1.85] 1.000 NA NA NA NA NA NA NA NA NA NA NA NA
Hu Daily persuasion utilized (partner’s view) 1.34*** [ 1.19, 1.54] 1.000 NA NA NA NA NA NA NA NA NA NA NA NA
Hu Daily pressure experienced 0.96 [ 0.70, 1.31] 0.595 NA NA NA NA NA NA NA NA NA NA NA NA
Hu Daily pressure utilized (partner’s view) 1.48* [ 1.04, 2.29] 0.987 NA NA NA NA NA NA NA NA NA NA NA NA
Hu Daily pushing experienced 0.95 [ 0.71, 1.31] 0.630 NA NA NA NA NA NA NA NA NA NA NA NA
Hu Daily pushing utilized (partner’s view) 1.32* [ 1.07, 1.68] 0.995 NA NA NA NA NA NA NA NA NA NA NA NA
Hu Day 0.86 [ 0.64, 1.15] 0.844 NA NA NA NA NA NA NA NA NA NA NA NA
Hu Actionplan 10.08*** [ 8.39, 12.18] 1.000 NA NA NA NA NA NA NA NA NA NA NA NA
Hu Daily weartime NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
Hurdle Between-Person Effects
Hu Mean persuasion experienced 1.39 [ 0.68, 2.81] 0.816 NA NA NA NA NA NA NA NA NA NA NA NA
Hu Mean persuasion utilized (partner’s view) 1.31 [ 0.64, 2.66] 0.770 NA NA NA NA NA NA NA NA NA NA NA NA
Hu Mean pressure experienced 0.41* [ 0.18, 0.92] 0.985 NA NA NA NA NA NA NA NA NA NA NA NA
Hu Mean pressure utilized (partner’s view) 0.50 [ 0.22, 1.14] 0.951 NA NA NA NA NA NA NA NA NA NA NA NA
Hu Mean pushing experienced 1.17 [ 0.41, 3.40] 0.619 NA NA NA NA NA NA NA NA NA NA NA NA
Hu Mean pushing utilized (partner’s view) 2.19 [ 0.76, 6.17] 0.928 NA NA NA NA NA NA NA NA NA NA NA NA
Hu Mean weartime NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
Random Effects
sd(Intercept) 0.31 [0.24, 0.41] NA 0.30 [0.23, 0.39] NA 0.60 [0.47, 0.78] NA 0.80 [0.47, 1.27] NA 1.16 [0.74, 1.74] NA
sd(Hurdle Intercept) 0.74 [0.57, 1.00] NA NA NA NA NA NA NA NA NA NA NA NA NA
sd(Daily persuasion experienced) 0.11 [0.08, 0.17] NA 0.05 [0.02, 0.08] NA 0.04 [0.00, 0.10] NA 0.17 [0.01, 0.43] NA 0.21 [0.01, 0.52] NA
sd(Daily persuasion utilized (partner’s view)) 0.08 [0.04, 0.13] NA 0.05 [0.03, 0.09] NA 0.07 [0.01, 0.13] NA 0.21 [0.01, 0.52] NA 0.48 [0.12, 0.98] NA
sd(Daily pressure experienced) 0.07 [0.00, 0.22] NA 0.04 [0.00, 0.13] NA 0.07 [0.00, 0.24] NA 0.55 [0.09, 1.16] NA 1.08 [0.13, 2.46] NA
sd(Daily pressure utilized (partner’s view)) 0.06 [0.00, 0.18] NA 0.03 [0.00, 0.11] NA 0.08 [0.00, 0.26] NA 0.42 [0.02, 1.58] NA 0.83 [0.04, 2.75] NA
sd(Daily pushing experienced) 0.09 [0.02, 0.17] NA 0.07 [0.01, 0.15] NA 0.05 [0.00, 0.14] NA 0.20 [0.01, 0.50] NA 0.25 [0.02, 0.61] NA
sd(Daily pushing utilized (partner’s view)) 0.08 [0.01, 0.15] NA 0.03 [0.00, 0.09] NA 0.07 [0.00, 0.17] NA 0.15 [0.01, 0.61] NA 0.26 [0.01, 0.96] NA
sd(Hu Daily persuasion experienced) 0.22 [0.04, 0.41] NA NA NA NA NA NA NA NA NA NA NA NA NA
sd(Hu Daily persuasion utilized (partner’s view)) 0.18 [0.02, 0.36] NA NA NA NA NA NA NA NA NA NA NA NA NA
sd(Hu Daily pressure experienced) 0.18 [0.01, 0.70] NA NA NA NA NA NA NA NA NA NA NA NA NA
sd(Hu Daily pressure utilized (partner’s view)) 0.23 [0.01, 0.90] NA NA NA NA NA NA NA NA NA NA NA NA NA
sd(Hu Daily pushing experienced) 0.58 [0.31, 0.99] NA NA NA NA NA NA NA NA NA NA NA NA NA
sd(Hu Daily pushing utilized (partner’s view)) 0.22 [0.01, 0.54] NA NA NA NA NA NA NA NA NA NA NA NA NA
Additional Parameters
sigma 0.68 [0.65, 0.70] NA 0.57 [0.56, 0.59] NA 0.96 [0.94, 0.98] NA NA NA NA NA NA NA 
report::report_system()

Analyses were conducted using the R Statistical language (version 4.4.1; R Core Team, 2024) on Windows 11 x64 (build 22635)

report::report_packages()
  • rmarkdown (version 2.29; Allaire J et al., 2024)
  • beepr (version 2.0; Bååth R, 2024)
  • R.methodsS3 (version 1.8.2; Bengtsson H, 2003)
  • R.oo (version 1.27.0; Bengtsson H, 2003)
  • R.utils (version 2.12.3; Bengtsson H, 2023)
  • brms (version 2.22.0; Bürkner P, 2017)
  • posterior (version 1.6.0; Bürkner P et al., 2024)
  • extrafont (version 0.19; Chang W, 2023)
  • digest (version 0.6.37; Eddelbuettel D, 2024)
  • Rcpp (version 1.0.13.1; Eddelbuettel D et al., 2024)
  • bayesplot (version 1.11.1; Gabry J, Mahr T, 2024)
  • lubridate (version 1.9.3; Grolemund G, Wickham H, 2011)
  • DHARMa (version 0.4.7; Hartig F, 2024)
  • fs (version 1.6.5; Hester J et al., 2024)
  • wbCorr (version 0.1.22; Küng P, 2023)
  • see (version 0.9.0; Lüdecke D et al., 2021)
  • tibble (version 3.2.1; Müller K, Wickham H, 2023)
  • patchwork (version 1.3.0; Pedersen T, 2024)
  • R (version 4.4.1; R Core Team, 2024)
  • openxlsx (version 4.2.7.1; Schauberger P, Walker A, 2024)
  • ggplot2 (version 3.5.1; Wickham H, 2016)
  • forcats (version 1.0.0; Wickham H, 2023)
  • stringr (version 1.5.1; Wickham H, 2023)
  • rvest (version 1.0.4; Wickham H, 2024)
  • tidyverse (version 2.0.0; Wickham H et al., 2019)
  • readxl (version 1.4.3; Wickham H, Bryan J, 2023)
  • dplyr (version 1.1.4; Wickham H et al., 2023)
  • purrr (version 1.0.2; Wickham H, Henry L, 2023)
  • readr (version 2.1.5; Wickham H et al., 2024)
  • xml2 (version 1.3.6; Wickham H et al., 2023)
  • tidyr (version 1.3.1; Wickham H et al., 2024)
  • ggridges (version 0.5.6; Wilke C, 2024)
  • knitr (version 1.49; Xie Y, 2024)
  • kableExtra (version 1.4.0; Zhu H, 2024)
report::cite_packages()